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Abstract — We consider distributed state estimation in a 
wireless sensor network without a fusion center. Each sensor 
performs a global estimation task — based on the past and current 
measurements of all sensors — using only local processing and 
local communications with its neighbors. In this estimation 
task, the joint (all-sensors) likelihood function (JLF) plays a 
central role as it epitomizes the measurements of all sensors. We 
propose a distributed method for computing, at each sensor, an 
approximation of the JLF by means of consensus algorithms. 
This "likelihood consensus" method is applicable if the local 
likelihood functions of the various sensors (viewed as conditional 
probability density functions of the local measurements) belong to 
the exponential family of distributions. We then use the likelihood 
consensus method to implement a distributed particle filter and 
a distributed Gaussian particle filter. Each sensor runs a local 
particle filter, or a local Gaussian particle filter, that computes a 
global state estimate. The weight update in each local (Gaussian) 
particle filter employs the JLF, which is obtained through the 
likelihood consensus scheme. For the distributed Gaussian parti- 
cle filter, the number of particles can be significantly reduced by 
means of an additional consensus scheme. Simulation results are 
presented to assess the performance of the proposed distributed 
particle filters for a multiple target tracking problem. 

Index Terms — Wireless sensor network, distributed state es- 
timation, sequential Bayesian estimation, consensus algorithm, 
distributed particle filter, distributed Gaussian particle filter, 
target tracking. 

I. Introduction 

Distributed estimation in wireless sensor networks has 
received significant attention recently (e.g., fl]-0). Appli- 
cations include machine and structural health monitoring, 
pollution source localization, habitat monitoring, and target 
tracking. Typically, a wireless sensor network is composed 
of battery-powered sensing/processing nodes — briefly called 
"sensors" hereafter — which possess limited sensing, computa- 
tion, and communication capabilities. 
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Centralized estimation techniques transmit sensor data to 
a possibly distant fusion center Q. This may require energy- 
intensive communications over large distances or complex 
multi-hop routing protocols, which results in poor scalability. 
Centralized techniques are also less robust, and less suitable 
if the estimation results have to be available at the sensors 
(e.g., in sensor-actuator networks [4|). Furthermore, the fusion 
center must be aware of the measurement models and, possi- 
bly, additional parameters of all sensors. By contrast, decen- 
tralized estimation techniques without a fusion center use in- 
network processing and neighbor-to-neighbor communications 
to achieve low energy consumption as well as high robustness 
and scalability. The sensors do not require knowledge of the 
network topology, and no routing protocols are needed. 

There are two basic categories of decentralized estima- 
tion techniques. In the first, information is transmitted in 
a sequential manner from sensor to sensor 0-G). In the 
second, each sensor diffuses its local information in an iterative 
process using broadcasts to a set of neighboring sensors (e.g., 
(8)). This second category is more robust but involves an 
increased communication overhead. It includes consensus- 
based estimation techniques, which use distributed algorithms 
for reaching a consensus (on a sum, average, maximum, etc.) 
in the network (9], [101 . Examples are gossip algorithms flTJI . 
consensus algorithms ATI , and combined approaches |fl2l . 

In this paper, we consider a decentralized wireless sensor 
network architecture without a fusion center and use consen- 
sus algorithms to perform a global estimation task through 
local processing and communications, in a way such that 
the final global estimate is available locally at each sensor. 
("Global" estimation means that the measurements of all 
sensors are processed by each sensor.) This can be based 
on the joint (all-sensors) likelihood function, abbreviated JLF, 
which epitomizes the measurements of all sensors. The JLF 
is then required to be known by all sensors. For example, a 
global particle filter (PF) |[T3l - lfT5l that processes all sensor 
measurements relies on the pointwise evaluation of the JLF to 
perform its weight update. 

The main contribution of this paper is a distributed method 
for calculating the JLF or an approximation of the JLF at 
each sensor. Generalizing our previous work in 1161 . 1171 . this 
method is suited to sensors with local likelihood functions 
that are members of the exponential family of distributions. 
A consensus algorithm — calculating sums — is used for a de- 
centralized, iterative computation of a sufficient statistic that 
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describes the (approximate) JLF as a function of the state 
to be estimated. Consequently, we refer to our method as 
likelihood consensus (LC). The LC scheme requires communi- 
cations only between neighboring sensors and operates without 
routing protocols. We furthermore propose an application 
of our LC method in a distributed PF scheme and in a 
distributed Gaussian PF scheme. Each sensor runs a local 
PF (or a local Gaussian PF lfl8l ) that computes a global 
state estimate incorporating all sensor measurements. At any 
given PF recursion, each local (Gaussian) PF draws a set of 
particles and updates their weights based on an evaluation of 
the JLF at these particles. For the distributed Gaussian PF, 
the number of particles employed by each local Gaussian PF 
can be significantly reduced by means of a second consensus 
scheme. 

Alternative consensus-based distributed PF schemes 
have been proposed in Ifl9l - ll24ll . The method described 
in lfl9l uses one consensus algorithm per particle to calculate 
products of local particle weights. To reduce the communi- 
cation requirements, the number of particles is kept small by 
an adaptation of the proposal distribution. Nevertheless, the 
number of consensus algorithms required can be significantly 
higher than in our approach. Furthermore, the random number 
generators of the individual sensors must be synchronized. On 
the other hand, since no approximation of the JLF is required, 
the performance can be closer to that of a centralized PF. The 
consensus-based distributed PFs proposed in l20l and ll2D rely 
on local PFs that update their weights using only the local 
likelihood functions instead of the JLF. Gaussian or Gaussian 
mixture approximations of local posteriors are then computed, 
and a consensus algorithm is used to fuse these approxi- 
mations. However, this fusion rule is suboptimal and leads 
to a performance loss. In J22), a novel gossiping approach 
implementing an approximation of the optimal fusion rule is 
employed to construct a Gaussian approximation of the global 
posterior. However, again only local likelihood functions are 
used by the local PFs, and the estimation performance is worse 
than in our approach. In ll23l . a distributed unscented PF is 
proposed that uses local measurements for proposal adaptation 
and an optimal consensus-based fusion rule to compute global 
estimates from local estimates. The distributed PF proposed in 
11241 operates across clusters of sensors and uses a modified 
maximum consensus algorithm to aggregate the local posterior 
distributions from all clusters. 

Distributed PFs that do not rely on consensus algorithms 
have been presented in fl25l - ll27l . In these methods, a path 
through the sensor network is adaptively determined by means 
of a decentralized sensor scheduling algorithm. Parametric 
representations of partial likelihood functions or of partial 
posteriors are transmitted along this path. The last sensor in the 
path obtains the complete global information and is thus able 
to compute a global estimate. In general, these methods are 
not as robust to sensor failure as the consensus-based methods. 
However, in certain applications, their communication require- 
ments may be much lower. 

This paper is organized as follows. In Section [II] we 
describe the system model and review sequential Bayesian 
estimation. To prepare the ground for the LC method, an 



approximation of the exponential class of distributions is 
discussed in Section [III] The LC method is presented in 
Section [IV] In Section [V] we consider the special case of 
additive Gaussian measurement noise. The application of LC 
to distributed particle filtering and distributed Gaussian particle 
filtering is considered in Section [VI] and IV11I respectively. 
Finally, in Section IVIIII the proposed distributed PFs are 
applied to multiple target tracking, and simulation results are 
presented. 

II. System Model and Sequential Bayesian 
Estimation 

We consider a wireless sensor network consisting of K 
sensors. At a given discrete time n, each sensor estimates a 
global M-dimensional state x n = (x n< \ ■ ■ ■ x„.m) T G K M 
based on all sensor measurements. The state evolves accord- 
ing to the state-transition probability density function (pdf) 
/(x„|x„_i). At time n, the fcth sensor (k £ {1,...,K}) 
acquires an N n> k -dimensional measurement z„.fc G E Ar " fc . The 
relationship between z„.fe and x„ is described by the local 
likelihood function^ f(z n ^\x n ), and the relationship between 
the all-sensors measurement vector z„ = (z„ ± - ■ ■ 2.^ K ) T and 
x„ is described by the JLF /(z„|x„). All z n> fc are assumed 
conditionally independent given x„, so that the JLF is the 
product of all local likelihood functions, i.e., 

K 

/(z„|x„) = Y[ f( z n,k\ x n)- (1) 
fe=l 

We write z± :n = (zj ■ ■ ■ z^) T for the vector of the measure- 
ments of all sensors up to time n. 

In the sequel, we will use the following assumptions. First, 
the current state x„ is conditionally independent of all past 
measurements, zi :n _i, given the previous state x„_i, i.e., 

/(x n |x„_i,zi : „_i) = /(x„|x„_i). (2) 

Second, the current measurement z„ is conditionally indepen- 
dent of all past measurements, S5i :n _i, given the current state 
x„, i.e., 

/(z„|x„,zi :n _i) = /(z„|x„). (3) 

Finally, sensor k knows the state-transition pdf /(x„|x„_i) 
and its own local likelihood function f(z n ^\x- n ) as well as 
the pdf /(xo) of the initial state xo, but it does not know the 
local likelihood functions of the other sensors, i.e., f(z nt k< |x„) 
for k'^k. 

We briefly review sequential Bayesian state estimation 
11281 . which will be considered as a motivating application of 
the LC method. At time n, each sensor estimates the current 
state x n from the measurements of all sensors up to time n, 
zi : „. For this task, we will use the minimum mean-square 
error (MMSE) estimator 1291 , 

^mmse a E { Xn | Zl:n } = j x „/(x„|z 1: „)dx„, (4) 

'The notation /(z n fc|x n ) suggests that x n is a random vector. However, 
for the LC method to be presented in Section IIVI x n is also allowed to 
be deterministic, in which case the notation /(z n fc;x n ) would be more 
appropriate. 
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which is implemented at each sensor. Here, a major problem — 
even in a centralized scenario — is to calculate the posterior 
pdf /(x„|zi : „). Using (0 and (0, the current posterior 
/(x„|zi : „) can be obtained sequentially from the previous 
posterior /(x n _i|zi : „_x) and the JLF /(z„|x„) by means of 
the following temporal recursion l28l : 



/(x„|Zl : „) 



/(z„|x„) / /(x n |x„_i)/(x„_i|zi : „_i)dx„_ 



/(Zn|zi:n-l) 



(5) 



However, for nonlinear/non-Gaussian cases, the computational 
complexity of sequential MMSE state estimation as given by 
© and (0 is typically prohibitive. A computationally feasible 
approximation is provided by the PF Ifl4l . fl31 . 1281 . In a PF, 

the (non-Gaussian) posterior /(x„ |zi „) is represented by a set 

(1) 

of samples (or particles) x„ , j = 1, . . . , J and corresponding 

( i) 

weights ra„ . 

As can be seen from (0 and (0, obtaining the global 
estimate x^f MSE at each sensor presupposes that each sensor 
knows the JLF /(z„|x„) as a function of the state x n (z„ 
is observed and thus fixed, and /(x n _i |zi :n _i) used in (0 
was calculated by each sensor at the previous time n — 1). 
In particular, a PF approximation of x?f MSE relies on the 

a) 

pointwise evaluation of the JLF at the particles x„ — i.e., 
on the evaluation of /(z rl |xl J ' ) ) — to obtain the weights Wn\ 
Since each sensor knows only its local likelihood function 
f(z n ,fc|x„), we need a distributed method for calculating the 
JLF at each sensor. Such a method is proposed in Section |IV] 
It is important to note that, although we consider distributed 
sequential Bayesian estimation and distributed particle filtering 
as a motivating application, the proposed method can also be 
used for other distributed statistical inference tasks that require 
the pointwise evaluation of the JLF at the individual sensors. 



III. Approximation of the Joint Likelihood 
Function 

The LC method can always be used if the local likelihood 
functions (viewed as conditional pdfs of the local measure- 
ments) belong to the exponential family of distributions. 
Typically, it requires an approximation of the local likelihood 
functions, and consequently of the JLF, which is discussed in 
the following. In Section IIV-BI we will consider a class of 
JLFs for which an approximation is not needed. 

A. Exponential Family 



In this paper, except in Section IIV-BI we assume that 
the local likelihood function of each sensor (viewed as the 
conditional pdf of z„ ^) belongs to the exponential family of 
distributions l30l. i.e., 



/(z„ jfe |x„) = Cn,k{i n ,k) exp(a 7lifc (x„)b„ ! fc(z„ ! fc) 

-d n ,fc(x„)), k = l,...,K, (6) 

with some time- and sensor-dependent functions c, l: k{-) £R+, 
a„, fc (-) £ R«, b„, fc (-) £ R«, and d„ jfe (-) £ R+, with arbitrary 
gSN. We furthermore assume that sensor k knows its own 
functions c n ^{-), a„, fc (-), b n , fe (-), and d n . k (-), but not c n ^(-), 



a„,fc'(-), b„, fc /(-), and d n , k >(-) for k' ^ k. Using (0, the JLF 
is obtained as 



A" 



/(z n |x„) = Y[ Cn,fe(zn,fe) exp (aj fc (x„) b„ ifc (z„ jfc ) 

k— 1 

- d„,fc(x n )) (7) 
= C„(z n )exp(5„(z n ,x n )) , (8) 

where 

K 

C„(z„) = JJc„ jfc (z„^) (9) 

and 



k=l 



K 

5„(z„,x„) = ^2 [ a X,fe( x n)b„,fe(z„, fc ) -rf„,fe(x n )] . (10) 
fc=i 

Note that the JLF (viewed as the conditional pdf of z„) 
also belongs to the exponential family. The normalization 
factor C„(z„) does not depend on the state x„ and is hence 
typically irrelevant; we will ignore it for now and consider 
it only at the end of Section IIV-AI Thus, according to (0, 
for global inference based on the all-sensors measurement 
vector z n , each sensor needs to know S n (z„ , x„) as a function 
of x„, for the observed (fixed) z„. However, calculation of 
5„(z„,x„) at a given sensor according to (TlOb presupposes 
that the sensor knows the measurements z„.fc and the functions 
a-n,fc( - )' b nj fe(-), and d n ^(-) of all sensors, i.e., for all k. 
Transmitting the necessary information from each sensor to 
each other sensor may be infeasible. 

B. Approximation of the Exponential Family 

A powerful approach to diffusing local information through 
a wireless sensor network is given by iterative consensus algo- 
rithms, which require only communications with neighboring 
sensors and are robust to failing communication links and 
changing network topologies 0111 . Unfortunately, a consensus- 
based distributed calculation of S n (z„,x n ) is not possible in 
general because the terms of the sum in dTOb depend on the 
unknown state x„. Therefore, we will use an approximation 
of S n (z n ,x. n ) that involves a set of coefficients not dependent 
on x n . This approximation is induced by the following ap- 
proximations of the functions a n fc(x„) and d n ,ki^n) in terms 



of given basis functions {<£ n ,r(xn)} r =i and {-0n,r(x„)} 
respectively: 



Rd 



Ra 



(Xri) ~ a n / e (x n ) ^ ^ Q^n.fc.r ^n.r (Xn ) 



(ID 



r=l 



= (x„) 



<4,fe(x„) = y^7n,fc,r^n,r(x w ) ■ (12) 
r=l 

Here, a rly k,r £ R 9 and jn^^ £ R are expansion coefficients 
that do not depend on x n . (For simplicity, the a n: k,r are 
referred to as coefficients, even though they are vector-valued.) 
The basis functions ^„ jr (x„) and ipn,r(x-n) do not depend on 
k, i.e., the same basis functions are used by all sensors. They 
are allowed to depend on n, even though time-independent 
basis functions may often be sufficient. We assume that sensor 
k knows the basis functions tp n r {-x. n ) and -0 njI .(x„), as well 
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as the coefficients ct n ,k,r and 7n,fc.r corresponding to its own 
functions a rai fe(x ra ) and d rai fc(x„), respectively; however, it 
does not know the coefficients of other sensors, a n! k', r and 
In.k'.r with k' 7^ k. The coefficients a n ,k,r and 7 n ,fe,r can 
either be precomputed, or each sensor can calculate them 
locally. A method for calculating these coefficients will be 
reviewed in Section UlI-CI 

Substituting a„. fc (x„) for a„, fc (x„) and J„. fc (x„) for 
d n ,k(x-n) m ( fTOt . we obtain the following approximation of 



Examples of basis functions tp n ,r(') and Vw(') are mono- 
mials (see the polynomial expansion discussed below), orthog- 
onal polynomials, and Fourier basis functions. The choice of 
the basis functions affects the accuracy, computational com- 
plexity, and communication requirements of the LC method. 

Example — polynomial approximation. A simple example 
of a basis expansion approximation ( fTTT i is given by the 
polynomial approximation 



M 



K 

E 

fe=i 

K 

E 

*:=i 



S„(z„,x„) = V] [aj jfc (x„)b„,fc(z„ ifc ) - J„, fc (x„)] (13) 



a„,fc 



(x„) = ^a„,fc,r J 



(17) 



r=l 

By changing the order of summation, we obtain further 

Ra Rd 

^i(Zn 7 X n ) — ^ -^-n^r (j^n) ^Pn^r (p^-ri) ~ T n _ r ^ n _ r (x n ) , 
r— 1 r— 1 

(14) 



with 



where r = (n- • • r^/) £ {0, . . . , i? p } M ; i? p is the degree of 
the multivariate vector-valued polynomial a„ fc(x„); ^ s 

short for ^=0" ' "X^Uo with the constraint J2m=i r ™ < 
R p ; and a n .k.r 6 is the coefficient vector associated with 
the basis function (monomial) V3„ jr (x rl ) = Ilm=i a; n,m (here, 
x n>m denotes the mth entry of x„). We can rewrite $1% in 
the form of (fTTT i by a suitable index mapping (ri ■ • • tm) £ 
{0,...,i? p } M o r G {l,...,i? Q }, where i? a = (\ +M ). 
An analogous polynomial basis expansion can be used" for 
dn,k(x-„) in ([T2l l. The polynomial basis expansion will be 
further considered in Section IV-BI 



K 



n.k) 



K 

P A C. Least Squares Approximation 

-L n,r — / J [n,k,r • 



k=l 



k=l 



(15) 

Finally, substituting 5 n (z n ,x n ) from (TBl i for 5„(z„,x ra ) in 
(O, an approximation of the JLF is obtained as 

/(z n |x„) oc exp(S„(z„,x„)) 



exp 



E 

r=l 



-Trijr 1pn t r (p^-n) 



(16) 



This shows that a sensor that knows A n r (z„) and r n r can 
evaluate an approximation of the JLF (up to a z n -dependent 
but x„-independent normalization factor) for all values of 
x„. In fact, the vector of all coefficients A„ iT .(z„) and r„ >r , 
t„(z„) = (4,i(zn) • • • A niRa (z n ) r n>1 • • • r„^ d ) , can be 
viewed as a sufficient statistic ll29l that epitomizes the total 
measurement z„ within the limits of our approximation. Be- 
cause of expression ( fT6l l, this sufficient statistic fully describes 
the approximate JLF /(z„|x„) as a function of x„. 

The expressions (fT4l and (fl3T l allow a distributed calcu- 
lation of 5„(z„,x„) and, in turn, of /(z„|x„) by means of 
consensus algorithms, due to the following key facts, (i) The 
coefficients A n r (z„) and r n r do not depend on the state 
x n but contain the information of all sensors (the sensor 
measurements z ra> fe and approximation coefficients a. n .k. r and 
ln,k,r for all k). (ii) The state x„ enters into S n (z n ,x n ) 
only via the functions if n ,r{-) and Vw(')> which are sensor- 
independent and known to each sensor, (iii) According to 1151 . 
the coefficients A n , r (z„) and r„ r are sums in which each term 
contains only local information of a single sensor. These facts 
form the basis of the LC method, which will be presented in 
Section HV^Al 



A convenient method for calculating the approximations 
a-n,fc(x n ) in STT[ and d^tfx,,) in (PL?! is given by least 
squares (LS) fitting 11311 - 1331 . We first discuss the calculation 
of the coefficients {ci!, l! A:,r}^ 1 of a ra> fc(x n ) at time n and 
sensor k. Consider J data pairs {(x^\,a n ^k(y^\))} J _ 1 , 

(i) ' ' 3 ~ 

where the state points x^ ; fe are chosen to "cover" those 
regions of the x n space where the JLF is expected to be 
evaluated when estimating x„. In particular, in the distributed 
PF application to be considered in Sections |VT] and IVHI the 
x « k rje tne predicted particles. With LS fitting, the 



coefficients ou 



,k,r 



are calculated such that the sum of the 



squared approximation errors at the state points xjri , i.e 



yj I 



Ja«,fc(x„; fe ) - a «,fc( x n,fe)ir' is minimized. 
To describe the solution to this minimization problem, we 

define the coefficient matrix Y n ^k = (at n ,k,i ■ ■ • ot n ,k,Ra) 

M. Ra * 9 , whose rows are the coefficient vectors {ot n ,k,r}^2i • 



r U) 



Furthermore, let 

V«4( x i'l) 

,m( x 2) 

A„, fc 4 (a n>k (x^ k ) ••• a„, fc (x^)) 



71, k 



T 



iJxq 



Then the LS solution for the coefficients {a n .fc.r} r=1 is given 
by ED 



Y„ 



(*n,fc*r 



Here, we assume that J>R a and that the columns of <f? ri! fc are 
linearly independent, so that ^Jfe^n.fc i s nonsingular. Note 
that J>R a means that the number of state points x^i is not 
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smaller than the number of basis functions <p n .r(x-n), for any 
given n and k. 

Similarly, the LS solution for the coefficients 
{in.kA^ti of 4i( x n) in <H3 is obtained as j nik = 
(*^ jfc *n.,fe) * T fe d nife , where 7„ >fe = (7„,fc,i • • • J n ,k,R d ) 
G R^, *„, fe G ' M Jxi ^ is defined like *„, fe but with 
{^n,r(-)}f=i replaced by {-0n,r(-)}f=r and d «,fc - 

{d n ,k( x n,k) "' d «,fc( x i J fc)) G K ' 7 - Here ' we assume that 
J>Rd and that the columns of ^ n ,k are linearly independent. 

To summarize, the number of state points must satisfy 

J > max{i? a , Rd] for any given n and fc. 

IV. Likelihood Consensus 

We now present the LC algorithm for local likelihood 
functions belonging to the exponential family, using the ap- 
proximation of the JLF discussed in Section Hill Subsequently, 
we will consider a class of JLFs for which an approximation 
is not needed. 

A. Distributed Calculation of the Approximate JLF - The LC 
Algorithm 

Based on the sum expressions (fL5T l. the sufficient statistic 
t n (z») = (Ai,i(z n ) • • • An >Ra (z n ) T„s ■ ■ ■ r n , Rd ) can be 
computed at each sensor by means of a distributed, iterative 
consensus algorithm that requires only communications be- 
tween neighboring sensors. Here, we use a linear consensus 
algorithm Hill for simplicity; however, other consensus algo- 
rithms (e.g., I34|) as well as gossip algorithms (e.g., ifTOl ) 
could be used as well. In what follows, the superscript W 
denotes the iteration index and A4 C {1, . . . , K}\{k} denotes 
a fixed set of sensors that are neighbors of sensor k. For 
simplicity, we only discuss the calculation of A n . r (z„), since 
the same principles apply to the calculation of r„ r in a 
straightforward manner. We explain the operations performed 
by a fixed sensor k; note that such operations are performed 
by all sensors simultaneously. 

At time n, to compute A, hr (z n ) = J2k'=x a Z,k',r 
x b n fe/(z n fc/), sensor k first initializes its local "state" as 
Cfe 0) — «I,fe jr bn,fc(Zn,fc)- This involves only the quantities 
z n ,fc> ctra,fc,r> and b Mj fc(-), all of which are available at sensor 
k; thus, no communication is required at this initialization 
stage. Then, at the ith iteration of the consensus algorithm 
(i G {1,2,...}), the following two steps are performed by 
sensor k: 

• Using the previous local state Cfc anc ' tne previous 
neighbor states , k' G A4 (which were received by 
sensor k at the previous iteration), the local state of sensor 
k is updated according to 

- w k,k^k + ^fe.fe'W 

The choice of the weights u) kk , is discussed in ll35l . Il36l . 
Here, we use the Metropolis weights 1361 

u<£ w =u k ,v = l + max{|A4|,|A/-H}' ^ ' 
[ 1 -T,k"eM k UJ k,k" , k'=k, 



where | Nk | denotes the number of neighbors of sensor k. 
(We note that knowledge at sensor k of |A4| and \Nk'\, 
k' G Nk is not required by certain other choices of the 
weights Il36l .) 

• The new local state (}^ is broadcast to all neighbors k' G 
Nk- 

These two steps are repeated in an iterative manner until a 
desired degree of convergence is reached. 

If the communication graph of the sensor network is 
connected, the state Q of each sensor k converges to the 
average J2k>=i ^k'^^k' (^n,k') = ~k A nA z n) as i^oo 
ifTTI . Therefore, after convergence, the states ^. i ~ >0 °' 1 of all 
sensors are equal and hence a consensus on the value of 
■h A nt r(z n ) is achieved. For a finite number i max of iterations, 

the states Ck""*^ W11 l ^ e (slightly) different for different sensors 
k and also from the desired value A n , r (z n ). In what 
follows, we assume that i max is sufficiently large so that 
_ A ntr (z n ) with sufficient accuracy, for all k. (In 
the simulations presented in Section [Villi i max G {7, 8, 9, 10}, 
which arguably does not imply impractical communication 
requirements.) Note that in order to calculate the coefficient 
A n ,r{z n ) from , each sensor needs to know K. This 

information may be provided to each sensor beforehand, or 
some distributed algorithm for counting the number of sensors 
may be employed (e.g., (37)). 

The consensus-based calculations of all A ntr (z n ), r = 
1, . . . , R a and all r„ >r , r = 1, . . . , Rd are executed simultane- 
ously, and their iterations are synchronized. These consensus 
algorithms taken together form the LC algorithm, which is 
stated in what follows. 



Algorithm 1: Likelihood Consensus (LC) 

At time n, the following steps are performed by sensor k (analogous 
steps are performed by all sensors simultaneously). 

1) Calculate the coefficients {a n ^,r}f2i anc ' {lri,k,r}f=i °f me 
approximations il lb and ( 112b . 

2) Consensus algorithm — A n , r (z n ): For each r = 1, . . . , R a '. 

a) Initialize the local state as (}^ = a^ fc r b nj fc(z„ ? fc). 

b) For i = 1,2, ...,i max (here, i max is a predetermined 
iteration count or determined by the condition that 
ICfe'-Cfc" 1 '! falls below a given threshold): 

(z) li) (i 1") 

• Update the local state according to Q — us), k Q 

I , ,(«) 

• Broadcast the new state Q to all neighbors k'EAfk- 

c) Calculate A n , r (z n ) = K( ( k tm) . 

3) Consensus algorithm — F n , r : For each r = 1, . . . , Rd- 

a) Initialize the local state as Q7 = j n ,k,r- 

b) Same as 2b). 

c) Calculate f„, r = KC£ W) . 

Finally, by substituting A nir (z n ) for A ntr (z n ) and r n . r for T n , r 
in J 1 6b . sensor k is able to obtain a consensus approximation of the 
approximate JLF /(z„|x n ) for any given value of x n . 
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Because one consensus algorithm has to be executed for 
each A 7hr (z n ), r = 1, . . . ,R a and r„ >r , r = 1, . . . , R d , the 
number of consensus algorithms that are executed simulta- 
neously is N c = R a + Rd- This is also the number of real 
numbers broadcast by each sensor in each iteration of the LC 
algorithm. It is important to note that R a and Rd do not depend 
on the dimensions N n ,k of the measurement vectors z„.fe, and 
thus the communication requirements of LC do not depend on 
the N, l: k- This is particularly advantageous in the case of high- 
dimensional measurements. However, R a and Rd usually grow 
with the dimension M of the state vector x„. In particular, if 
the MD basis {%,rW} fl =i is constructed as the M-fold 
tensor product of a ID basis {<f n ,r(x)}?l v then R a = R^f , 
and similarly for the MD basis {VVi,r(x n )} _f r 

So far, we have disregarded the normalization factor 
C n (z n ) occurring in ([8). If this factor is required at each 
sensor, it can also be computed by a consensus algorithm. 
From (|9), 



logC„(z n ) 



K 

E 

fc=i 



log c n , fe (z 



Since this is a sum and c n ^k{z n ,k) is known to each sensor, 
a consensus algorithm can again be used for a distributed 
calculation of logC n (z„). 

B. Distributed Calculation of the Exact JLF 

The basis expansion approximations (fTTT > and (fT2l can 
be avoided if the JLF /(z„|x„) has a special structure. In 
that case, the exact JLF can be computed in a distributed 
way, up to errors that are only due to the limited number 
of consensus iterations performed. We note that the special 
structure considered now is compatible with the exponential 
family structure considered so far, but it does not presuppose 
that structure. 

Let t„(z„) = (t„,i(z„) • • • t„,p(z„)) be a sufficient 
statistic for the estimation problem corresponding to /(z„ |x„). 
According to the Neyman-Fisher factorization theorem 
/(z„|x„) can then be written as 



/(z„|x n ) = /i(z n )/ 2 (t n (z n ),x n ) 



(18) 



Typically, the factor /i(z ra ) can be disregarded since it does 
not depend on x„. Thus, t„(z„) epitomizes the total mea- 
surement z n , in that a sensor that knows t„(z„) and fi{- , •) 
is able to evaluate the JLF f(z n |x„) (up to an irrelevant factor) 
for any value of x n . Suppose further that the components of 
t„(z„) have the form 



K 



k=l 



r ]n,k,p{Z'n,k) 



P 



1, 



,P, 



(19) 



with arbitrary functions r) n! k, p (-), and that sensor k knows its 
own functions r) n! k, p (-) but not i] n .k'.p{-), k'^k. Based on the 
sum expression ( fT9l . we can then use consensus algorithms 
as described in Section IIV-AI with obvious modifications, to 
calculate t„(z„) and, thus, the JLF /(z„|x n ) in a distributed 
manner. 

Clearly, an example where exact calculation of the JLF is 
possible is the case where /(z„|x„) belongs to the exponential 



family (0, with functions a ni fc(x„) and d nj fc(x„) that can 
be exactly represented using expansions of the form (fTTT > and 
(H3, i.e., a„ jfc (x n ) = J2r=i a Jl ,fc,r<^n,r(x„) and d„,fc(x„) = 
12r=i 7«,fc,r'0n,r(x n ). This is a special case of ( TT8l and ( fT9l . 
with (cf. 

/ 2 (t„(z 



where P — R a + Rd and 



i Rat 



- n.p—R a 



E 



K 

k=l ln,k,p-R a i 

p = Ra + l, 



Pn,p (p^-n ) 



<Pn,p(x„), p=l,...,R a , 
v '0n,p-fl a (x n ), p = R a + l,...,P. 

Equivalently, i„ iP (z n ) is of the form ( fT9l . with (cf. ( fT5] >) 

a Ifc,p b ",fc(z„. fc ), p=l,...,i2 a , 



Vn,k,p{ z n,k) 



' ^fn^k^p—Ra i 



P = R a +l, 



,P. 



V. Special Case: Gaussian Measurement Noise 

In this section, we consider the important special case of 
(generally nonlinear) measurement functions and independent 
additive Gaussian measurement noises at the various sensors. 
We will also develop the application of the polynomial ap- 
proximation that was briefly introduced in Section IIII-BI 

A. Measurement Model 

The dependence of the sensor measurements z ra> fc on 
the state x„ is described by the local likelihood functions 
/(z„ fc|x„). Let us now assume, more specifically, that the 
measurements are modeled as 

z«,fc = h n> fc(x n ) + v„ jfc , k = 1, . . . , K , (20) 

where h nJ t(-) is the measurement function of sensor k and 
v n ,fc ~ jV(0, Q n ,fc) is zero-mean Gaussian measurement noise 
that is independent of x„/ for all n'. We furthermore assume 
that v„.fc and v n /fe/ are independent unless (n, k) = (n', k'). 
Under these assumptions, the z n> k are conditionally indepen- 
dent given x n , i.e., ([TJi holds. The local likelihood function of 
sensor k is here given by 

/(z njfc |x n ) 

= c^fcexpf- i[z n;fc -h„ ;fc (x n )] T Q n ;j £ [z nifc -h nife (x n )]^ , 

(21) 

with c ni k — [{2Tt) Nn - k det{Q ri! fc}]~ 1//2 . Furthermore, using 
([TJ, the JLF is obtained as 




1 K 



z„.fe-h„. fe (x„)] T Q n ^. [z„,i 



k=l 



•h„ lfe (x„)] 
(22) 
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with On = Cn,k- the directly and indirectly obtained coefficients 7 ra> fc jr will be 

The local likelihood function /(z„ fc|x„) in (f2TT> is a different. Furthermore, if the indirectly obtained coefficients 
special case of the exponential family (O, with are used, the approximate JLF /(z n |x ra ) is a valid pdf in 



a„,fc(x„) = h„ ifc (x„), (23) 



the sense that J f(z n \x n )dz n = 1 holds exactly, not only 
approximately. The number of consensus algorithms that are 
b n ,k(zn,k) = Q n ^ z n,k , executed is N c = R a + Rd = R a + R\. Again, this does not 

/ 1 T , \ depend on the dimensions N n> k of the measurement vectors 

c«.,fc( z n,fc) = Cn,k cxp I - - z nk Q nk z„,k I , Znk smce ^ a joes not depend on N njk . 

<fc(x„) = -hJ ife (x n )Q- 1 fe h n)fc (x n ). (24) 



2 ' B. Polynomial Approximation 

Consequently (see (TTOli). 

if 

S n (z n ,X n ) = E h^.fc ( X » ) Qn.l- 



fc=l 



The polynomial approximation was introduced in Section 
IIII-BI We will now apply it to the case of Gaussian measure- 
(25) ment noise studied above. Using ( fTTI i. we obtain for 



We now approximate a„. fe (x„) and d n , fc (x n ) by truncated - , , _ r , * _ \p rr Tm nQ , 
basis expansions a nj fc(x n J and d„ i fe(x n J of the form (II 11 1 and 
(fl2l) . respectively. According to (l23l . approximating a n ,fc(x n ) 

is equivalent to approximating the sensor measurement func- Inserting this into (f27]i yields 

tion h„.fc(x„) (which is also the mean of /(z ni fc|x„) in (fJTJ). 2 _r m 

ThUS ' fla < fc (x„) = J2ln,k,rf[x r n7m , (30) 



r— m— 1 



with 



(31) 



an,fe(x n ) = h n! fc(x„) = ^2 a n,k,T <£n,r(x n ) ■ (26) 
r=l 

Furthermore, an approximation of e?„.fc(x„) of the form (fT2l 1 j j 

can be obtained in an indirect way by substituting in ( |24T i the 2 ^-^ n, * ! ' r ' ™> fc 1 

above approximation h„j £ (x TI ) for h„ fc(x„); this yields r'+r" = r 

j /„ \ _ 1 CT / \ n -i r / \ / 97 x Next, inserting expressions d29| i and (l30l into (TTBl , we obtain 

"n.fcl.Xn; — - n nk [x. n ) \4 nk n n ,k\ x -n) 

Ra Ra K 2Rp M 

fe=l r=0 m=l 



2 

ri=l r 2 = l 



(28) with 



Using a suitable index mapping (n,r 2 ) G {l,...,i? a }x J aj fe r b„ !fc (z„. fc ) - 7 n ,fc,r, reT^i 

{l,...,Ra} O r G {l,...,i?4, we can write (f28j in the *Wl*n,*J ~ i_ 7nfcr , re ft 2 , 

form G2) 



where TZi is the set of all r = (n- • • tm) G {0, . . . , i? p } M 
d„,fc(x„) = 2^7„, fc ,r'0r l ,r(x„) , such ^ £^" =1 r m < i? p and ft 2 is the set of all r 6 

r=1 {0, . . . , 2i? p } M \ Hi such that £^ =1 r m < 2i V Finally, 

with i? d = R 2 a , j n ,k, r = 2 £ *n,ife,riQn,fc a n>fe>»'a' anc ^ Vwfxn) changing the order of summation in d32l gives 
= (pn.r 1 (xn) ( Pn,r 2 ( :>c n)- It is easily verified that with this spe- 



2-Rp M 



cial basis expansion approximation of d n ,fe(x„), the resulting ~ d \ tt r m 

approximate JLF can be written as Mz n , x «) = 2^ B nA z n) [[ x™ m , (34) 



r— rn— 1 



/(z n |x„) with 



C„exp( -- ^[z„ ifc -h„ jfe (x„)] T Q n 1 fc [z„,fc-h„ j fe(x n )] J , B„, r (z„) = ^j9 w ,fc,r(Zn,A 

V fc=i / fc=i 



(35) 



which is d22l i with h n ^(x„) replaced by h„^(x n ). This It should be noted that (l34l is a special case of ( fT4l i. The 

means that only the mean of /(z„|x„) is changed by this coefficients J3„ jr (z„) can again be calculated using a consen- 

approximation. sus algorithm. For each time n, the number of coefficients 

In the additive Gaussian noise setting considered, the LC B nyT (z n ), and hence the number of consensus algorithms that 

method operates almost as in the general case. The only have to be executed in parallel, is given by N c = ( 2i ^ M ) — 1. 

difference is in Step 1 of Algorithm 1: instead of calculating Here, the subtraction of 1 is due to the fact that the coeffi- 

the coefficients 7„jt. r directly, using, e.g., a separate LS fitting, cient B n T= o(z n ) need not be calculated: according to ( f34b . 

we obtain them in an indirect way as described above. Hence, B„,o(z n ) corresponds to a JLF factor that does not depend on 

the computational complexity is reduced. Note that in general, x n and is hence irrelevant. 
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VI. Distributed Particle Filtering 

In this section, we show how the LC method can be applied 
to obtain a distributed PR By way of preparation, we first 
review a standard centralized PF lfl4l . 03], J28). 

A. Review of Centralized Particle Filtering 

The centralized PF is implemented at a fusion center 
that knows the all-sensors measurement vector z n and the 
functional form of the JLF /(z n |x„). The PF maintains a set 
of samples (or particles) {x^ and associated weights 

{wn '}^, which establish the following approximative sam- 
ple representation of the posterior pdf /(x„|zi : „): 

J 

/ 5 (x„|zi :n ) = ^W^^Xn-xW). 

3=1 

The MMSE estimate in (|4]i can then be approximated by the 
mean of /,5(x„|zi :n ), which is equivalent to a weighted sample 
mean: 

./ 



E 

3 = 1 



wl j) x« 



(36) 



At each time step n, when the new measurement vector z ra 
becomes available, new particles and weights are calculated 
by a PF algorithm that is based on the recursion (0. 

Many PF algorithms have been proposed lfT3l - lfT31 . l28l . 
Here, we consider a sequential importance resampling filter 

IfTJI . ITT51 . which performs the following steps. For initial- 

(i) 

ization (n = 0), J particles Xq are sampled from a prior 
distribution /(xo), and the weights are set to Wq^ = l/J. 
Then, three steps — resampling, sampling, and weight update — 
are repeated for every n. In the resampling step, J resampled 
particles x^2i are obtained by sampling with replacement 
from the set of previous particles {x^_\ } J -,_ x , where the prob- 
ability of sampling x^ ^ is w^\. In the sampling step, for 
each resampled particle x„_ 1; a new, "predicted" particle x„ 
is sampled from ,/(x Il |x^'l 1 ), i.e., from the state-transition pdf 
/(x„|x„_i) evaluated at x„_i = x^ij^. In the weight update 
step, the weight associated with each particle is calculated 
as 



w 



(j) 



r 0')^ 



ELi/(^|xf } ) 



(37) 



This requires each sensor to know the JLF /(z„|x n ) as a 
function of the state x„, because the weight update in (|37| | 
requires the pointwise evaluation of the JLF. Therefore, an 
approximation of the JLF is provided to each sensor in a 
distributed way by means of the LC method. No routing of 
measurements or other sensor-local data is needed; each sensor 
merely broadcasts information to neighboring sensors. The 
algorithm is stated as follows. 

Algorithm 2: LC-based Distributed PF (LC-DPF) 

At time n, the local PF at sensor k performs the following steps, 
which are identical for all k. (Note that these steps are essentially 
analogous to those of the centralized PF of Section fVI-AI except that 
an approximation of the JLF is used.) 

1) At the previous time n- 



x i-i.fc ^ weights w^ ] _ lk 



1, sensor k calculated J particles 
which together represent the pre- 
vious global posterior /(x„_i |zi :n _i). The first step at time n 
is a resampling of { (xi 
resampled particles x 



0) 

71- 

(j) 



,,(1)^ fc ) } J _ , which produces J 
Here, the x"_ 1 k are obtained by 



sampling with replacement from the set {x^j^ k } J ,_ 1 , where 
x n-i k ls sam pled with probability w^_\ k . 

2) For each x^ij k , a new, "predicted" particle x^' fe is sampled 
from /(x n |x n _i)j _ (j) . 

3) An approximation /(z„ |x n ) of the JLF /(z n |x n ) is computed 
by means of LC as described in Section IIV-AI This step 
requires communications with neighboring sensors. The local 
approximation at sensor k can be calculated by means of 
LS fitting as described in Section IIII-CI using the predicted 



particles {^ k } J 



4) The weights associated with the predicted particles x ?l fe 
tained in Step 2 are calculated according to 



ob- 



„0) 



E^ =1 /WX^ 



3 = 1,..., J. 



(38) 



This involves the approximate JLF /(z n |x n ) calculated in Step 
3, which is evaluated at all predicted particles x! 



5) From { (x. 



CO . 



,0) 



)} J _ 1 , an approximation of the global 



MMSE state estimate $4$ is computed according to yb), i.e.. 



x„ t 



= E 



w ( JWJ\, 



The recursion defined by Steps 1-5 is initialized at n = by 
J particles Xq 3 k sampled (at each sensor) from a suitable prior pdf 

CO 



Finally, the state estimate x„ is calculated from {(x?, /(x ), and by equal weights w ' k = l/J 



w ^)}'j=i acc °rding to d36l l. 

B. Distributed Particle Filtering Using LC 

Next, we develop a distributed implementation of the 
sequential importance resampling filter reviewed above, in 
which each sensor acts similarly to the fusion center of the 
centralized PF. More specifically, sensor k tracks a particle 
representation of the global posterior /(x„|zi : „) using a local 
PF. For each n, it obtains a state estimate x„jt that is based 
on zi : „, i.e., the past and current measurements of all sensors. 



Through the above recursion, each sensor obtains a global 
quasi-MMSE state estimate that involves the past and current 
measurements of all sensors. Because of the use of LC, this 
is achieved without communicating between distant sensors or 
employing complex routing protocols. Also, no particles, local 
state estimates, or measurements need to be communicated 
between sensors. The local PF algorithms running at different 
sensors are identical. Therefore, any differences between the 
state estimates x n fe at different sensors k are only due to 
the random sampling of the particles (using nonsynchronized 



9 



local random generators) and errors caused by insufficiently 
converged consensus algorithms. 

C. Communication Requirements 

We now discuss the communication requirements of our 
LC-based distributed PF (LC-DPF). For comparison, we also 
consider the centralized PF (CPF) of Section IVI-AI in which 
all sensor measurements are transmitted to a fusion center 
(FC), and a straightforward distributed PF implementation (S- 
DPF) in which the measurements of each sensor are transmit- 
ted to all other sensors. Note that with the S-DPF, each sensor 
performs exactly the same PF operations as does the FC in 
the CPF scheme. 

For the CPF, communicating all sensor measurements at 
time n to the FC requires the transmission of a total of 
J2k=i HkN nt k real numbers within the sensor network l25l . 
Here, Hk denotes the number of communication hops from 
sensor k to the FC, and N n; k is the dimension of z„fc. 
Additional information needs to be transmitted to the FC if 
the FC does not possess prior knowledge of the JLF Finally, 
if the state estimate calculated at the FC is required to be 
available at the sensors, additional MH' real numbers need to 
be transmitted at each time n. Here, H' denotes the number of 
communication hops needed to disseminate the state estimate 
throughout the network. A problem of the CPF using multihop 
transmission is that all data pass through a small subset of 
sensors surrounding the FC, which can lead to fast depletion 
of the batteries of these sensors. 

With the S-DPF, disseminating the measurements of all 
sensors at time n to all other sensors requires the transmission 
of J2k=i H'^N n k real numbers 11251 . where H' k ' is the number 
of communication hops required to disseminate the measure- 
ment of sensor k throughout the network. Again, additional 
information needs to be transmitted if the JLF is not known 
to all sensors. 

Finally, the proposed LC-DPF requires the transmission of 
KIN C real numbers at each time n, where / is the number of 
consensus iterations performed by each consensus algorithm 
and N c = R a + Rd is the number of consensus algorithms 
executed in parallel (see Section ITV- At . In contrast to the CPF 
and S-DPF, this number of transmissions does not depend 
on the measurement dimensions N n y, this makes the LC- 
DPF particularly attractive in the case of high-dimensional 
measurements. Another advantage of the LC-DPF is that 
no additional communications are needed (e.g., to transmit 
local likelihood functions between sensors). Furthermore, the 
LC-DPF does not require multihop transmissions or routing 
protocols since each sensor simply broadcasts information to 
its neighbors. This makes the LC-DPF particularly suited to 
wireless sensor networks with dynamic network topologies 
(e.g., moving sensors or a time-varying number of active 
sensors): in contrast to the CPF and S-DPF, there is no need to 
rebuild routing tables each time the network topology changes. 

On the other hand, the computational complexity of the 
LC-DPF is higher than that of the S-DPF because the approx- 
imation described in Section UTT1 needs to be computed at each 
sensor. Overall, the LC-DPF performs more local computa- 
tions than the S-DPF in order to reduce communications; this 



is especially true for high-dimensional measurements and/or 
high-dimensional parametrizations of the local likelihood func- 
tions. Since the energy consumption of local computations is 
typically much smaller than that of communication, the total 
energy consumption is reduced and thus the operation lifetime 
is extended. This advantage of the LC-DPF comes at the cost 
of a certain performance loss (compared to the CPF or S-DPF) 
due to the approximate JLF used by the local PFs. This will 
be analyzed experimentally in Section [Villi 

VII. Distributed Gaussian Particle Filtering 

Next, we propose two distributed versions of the Gaussian 
PF (GPF). The GPF was introduced in ffT~8l as a simplified 
version of the PF using a Gaussian approximation of the pos- 
terior /(x ra |zi : „). The mean and covariance of this Gaussian 
approximation are derived from a weighted particle set. The 
particles and their weights are computed in a similar way as 
described in Section [VI] with the difference that no resampling 
is required. This results in a reduced complexity and allows 
for a parallel implementation ll38l . 

A. Distributed Gaussian Particle Filtering Using LC 

In the proposed distributed GPF schemes, sensor k uses 
a local GPF to track the mean vector /i.„ & and covariance 
matrix G n> k of a local Gaussian approximation Af(fi n ,k, C n ,k) 
of the global posterior /(x„|zi : „). The state estimate x„ ^ of 
sensor k at time n is defined to be the current mean, i.e, 
x n> fe = Hn,k- The calculation of this estimate is based on the 
past and current measurements of all sensors, zi :n . As with 
the distributed PF described in Section IVI-BI (Algorithm 2), 
these measurements are epitomized by an approximation of 
the JLF, which is provided to each sensor by means of LC. A 
statement of the algorithm follows. 



Algorithm 3: LC-based Distributed GPF (LC-DGPF) 

At time n, the local GPF at sensor k performs the following steps, 
which are identical for all k. 

1) J particles {x^i k}j-i are sam pl e d from the previous local 
Gaussian approximation M{n n -i,k, C„_i,fc), where (J, n -i,k 
and C n -i,k were calculated at time n — 1. Note that this 
sampling step replaces the resampling step of the distributed 
PF of Section 171-51 (Step 1 in Algorithm 2). 

2)— 4) These steps are identical to the corresponding steps of the 
distributed PF of Section IVI-BI (Algorithm 2); they involve 
LC (Step 3) and result in a set of "predicted" particles and 
corresponding weights, {{^ k ' w n,k)}j=i- 

5) From { (x„ k , w!^\) tne mean Mn.fc and covariance C n ,fc 
of the Gaussian approximation M(fJ, n ,k,C nt k) of the current 
posterior /(x n |zi ;n ) are calculated as 

j 

Mn.fc = J2 W n?k*n?k 
;' = 1 

J Q9) 
r , _ NP,„«) „W)„0)T_ T 
^n,k — 2^ W n,fc X n,fe X n,fc Mn,t ■ 

J = l 
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The state estimate x„,k (approximating xJf MSE in l[4]l) is taken 
to be the posterior mean /Ltn.fe. 
The recursion defined by Steps 1-5 is initialized as in Algorithm 

2. 



B. Reduced-Complexity Method 

We next present a reduced-complexity variant of the LC- 
DGPF described above, in which each of the K local GPFs 
uses only J' = J/K particles. Here, J is chosen such that J' 
is an integer and J' > max{_R a , Rd] (cf. Section ITlI-Ct . The 
sets of J' particles of all local GPFs are effectively combined 
via a second stage of consensus algorithms, such that a "virtual 
global GPF" with J = KJ' particles is obtained. In other 
words, J particles — which, in the LC-DGPF, were used by 
each individual sensor separately — are "distributed" over the 
K sensors. As a consequence, the computational complexity 
of the local GPFs is substantially reduced while the estima- 
tion accuracy remains effectively unchanged. This advantage 
comes at the cost of some increase in local communications 
due to the additional consensus algorithms. 

This reduced-complexity method is similar to a parallel 
GPF implementation proposed in Il38l . which uses multiple 
processing units — corresponding to our sensors — collocated 
with a central unit. However, instead of a central unit, we 
employ distributed consensus algorithms to combine the partial 
estimates (means) and partial covariances calculated at the 
individual sensors. Another difference from l38l is the use 
of an approximate JLF that is obtained in a distributed way 
by means of LC. The algorithm is stated as follows. 



Algorithm 4: Reduced-complexity LC-DGPF 
(R-LC-DGPF) 

At time n, the local GPF at sensor k first performs Steps 1-3 of 
the LC-DGPF algorithm described in Section IVII-AI (Algorithm 3), 
however using J' = J/K rather than J particles. The remaining 
steps, described in the following, are modified versions of Steps 4 
and 5 of Algorithm 3, as well as an additional consensus step. 

4) Nonnormalized weights are calculated as (cf. ((38} ) 

<l = /Wx^), 3 = 1,..., J'. 

This requires evaluation of the approximate JLF /(z„|x„), 
which was calculated in Step 3 using LC, at the J' predicted 
particles {x„ drawn in Step 2. Furthermore, the sum of 

the J' nonnormalized weights is computed: 

j' 

3=1 

5) From the weighted particles { (x„ fc , w^\) a partial 

nonnormalized mean and a partial nonnormalized correlation 

are calculated as 

j' J' 
.,' _V^,r,0')„W) r>' _ \p,r,0) v 0) Y (i)T 

Mn.fc - 2^^n,fc X n,fe ' - / , W n.k X ».fc X «,fc > 

j=l j=l 

(40) 

respectively. Note that Steps 4 and 5 are carried out locally at 
sensor k. 



6) The partial means and correlations from all sensors are com- 
bined to obtain the global mean and covariance: 

1 r 1 A , 

Wn k=l W " k=l 

(41) 

where 

W n = Y^W n ,k (42) 

fc=i 

is the global sum of all particle weights. The sums over all 
sensors in J41t and ( 1421 ) are computed in a distributed manner 
by means of consensus algorithms. The normalization by W n 
and subtraction of UnfJ-Z in d41b are performed locally at each 
sensor after convergence of these consensus algorithms. The 
state estimate x n is taken to be fi n . 



As a result of this algorithm, all sensors obtain identical 
x n = A*n and C„ provided that the consensus algorithms 
are sufficiently converged. Therefore, we omit the subscript 
k indicating the sensor dependence (cf. fl39l)), i.e., we write 
x n = A 4 " instead of x„,fc = fi n ,k an d C„ instead of C n ,h for 
all k. 

It is easily seen from (l40b — (|42T > that fi n and C n are actually 
the result of an averaging (summation) over J particles (note 
that J' = J/K particles are sampled independently at each of 
the K sensors). Therefore, under the assumption that the con- 
sensus algorithms used to calculate the sums over all sensors 
in (f4TT > and d42l i are converged, fi n and C„ should ideally be 
effectively equal to the corresponding quantities obtained by 
the LC-DGPF. However, a certain performance degradation 
is caused by the fact that the LS fitting performed at each 
sensor (see Section Ull-Cb is now based on only J' = J/K 
predicted particles x^ fc , and hence the resulting approximate 
local likelihood functions and, in turn, the approximate JLF 
will be less accurate. In Section [Villi we will show by means 
of simulations that this degradation is very small. 

C. Computational Complexity and Communication Require- 
ments 

We compare the computational complexity and commu- 
nication requirements of the LC-DGPF and of its reduced- 
complexity variant discussed above (abbreviated R-LC- 
DGPF). We will disregard Steps 2 and 3 of the LC component 
(Algorithm 1), because their complexity and communication 
requirements are identical for the LC-DGPF and R-LC-DGPF; 
furthermore, their complexity is typicall}@ much lower than 
that of the remaining steps (local GPF algorithm and LS 
approximation). 

The complexity of the local GPF algorithm and of the 
LS approximation in the LC scheme (Step 1 of Algorithm 
1) depends linearly on the number of particles OTI . Il38ll . 
Thus, reducing the number of particles at each sensor from 
J to J' = J/K reduces this complexity by a factor of K. It 
follows that the R-LC-DGPF is significantly less complex than 

2 The complexity of Steps 2 and 3 of Algorithm 1 is linear in the number 
of consensus algorithms and in the number of consensus iterations; these 
numbers depend on the specific application and setting. 
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the LC-DGPF. (The complexity of the additional consensus 
algorithms required by the R-LC-DGPF is typically negligible 
compared to the other operations.) The additional communica- 
tion requirements of the R-LC-DGPF relative to the LC-DGPF 
are determined primarily by the speed of convergence (i.e., 
number of iterations /) of the additional consensus algorithms, 
which depends mainly on the second smallest eigenvalue of 
the Laplacian of the communication graph ll39l . and by the 
state dimension M. More specifically, the additional number 
of real numbers transmitted in the entire sensor network at 
each time n is KIN' C , where N' c = M + M(M+ 1)/2 + 1 is 
the number of additional consensus algorithms, i.e., of (scalar) 
consensus algorithms needed to calculate the mean vector and 
covariance matrix in (f4lT> as well as the total weight in d42l) . 
Since N' c is of order M 2 , the R-LC-DGPF has a disadvantage 
for high-dimensional states. 

The reduced operation count of the R-LC-DGPF relative to 
the LC-DGPF can be exploited in two alternative ways, which 
represent a tradeoff between latency and power consumption. 
First, the processing time can be reduced; this results in a 
smaller latency of the R-LC-DGPF relative to the LC-DGPF, 
provided that the delays caused by the additional consensus 
algorithms are not too large. Thus, the R-LC-DGPF may be 
more suitable for real-time applications; however, the power 
consumption is higher due to the increased communications. 
Alternatively, if latency is not an issue, the processor's clock 
frequency can be reduced. The processing time can then be 
made equal to that of the LC-DGPF, while the processor's 
power consumption is reduced due to the lower clock fre- 
quency l40l . Thereby, the overall power consumption of the 
R-LC-DGPF is smaller relative to the LC-DGPF, provided 
that the additional power consumption due to the increased 
communications is not too large. However, the total latency 
is increased by the delays caused by the additional consensus 
algorithms. 

VIII. Numerical Study 

We will now apply the proposed LC -based distributed PF 
algorithms to the problem of tracking multiple targets using 
acoustic amplitude sensors. We will compare the performance 
of our methods with that of the centralized PF and state-of- 
the-art distributed PFs. 

A. Acoustic-Amplitude-Based Multiple Target Tracking 

We consider P targets (P assumed known) moving inde- 
pendently in the x-y plane. The pth target, p E {1, . . . , P}, is 
represented by the state vector = (x^ ±i P ' Vn^) 
containing the target's 2D position and 2D velocity. The 
overall state vector is defined as x„ = (x„ • • • xi P ' T ) . 



Each vector x„ evolves independently of the other vec 

-GO 
p-^n—l 



; x£f ' according to x„ = G„xf p ^ 



W p ui p) . I 



~ Af(Q2, is Gaussian driving noise, with u„ and 



,0) 

u^ p ^ independent unless (n,p) = (n',p'), and G p 6 W lxi 
and W p G R 4x2 are system matrices that will be specified 
in Section IVIII-BI This model is commonly used in target 
tracking applications JTSJ, BT1 - B31 . It follows that the overall 



state vector x„ evolves according to 
x n = Gx„_i + Wu„ , 



1,2, 



where G = diag{Gi, . . . , G P }, W = diag{Wi, . . . , W P }, 
and u n 4 (ul 1)T - • • ul P)T ) T ~ AT(0 2P , a 2 u I 2P ). 

Each target p emits a sound with a (root mean-squared) 
amplitude A p that is assumed constant and known. At the 
position of sensor k, denoted £ n ,k, the sound amplitude due 
to target p is modeled as A p /\\pn — £ n> k\\ K , where pi P ' = 
(xn Vn) is the position of target p and k is the path 
loss exponent BTIl . fi4ll . B31l . The (scalar) measurement z n .t 
obtained by sensor k at time n is then given by 



Zn.k — ^ni(^n) "F ^n,k ; 

P 

with h nj k(Xn) = ^2 



(43) 



where v n ,k ~ jV(0, a 2 ) are zero-mean Gaussian measurement 
noise variables of equal variance a 2 . We assume that v n> k is 
independent of x n / for all n', and that v n> k and u n ',fc' are inde- 
pendent unless (n, k) = (n', k'). Note that this measurement 
model is a special instance of (l20l i, and that z rlJ t does not 
depend on the velocities i„ and y„ . The local likelihood 
functions and the JLF are respectively given by (cf. (I2TI 1. (l22l ) 



/(^n.fclXn) 



/(Zn|x„) 



exp 



2o* 



[z„,fc-/in,fc(x„)] 2 (44) 



and hence (cf. (T25] l) 



exp 



1 K 



= (x„)] 5 



k=l 



1 K 

v k=i 



with /i„,fc(x„) given by (|43j. 

In general, the sensor positions are allowed to change 
with time n. (However, we used static sensors for simplicity.) 
Each sensor is supposed to know its own position but not the 
positions of the other sensors. The sensor positions (which are 
contained in the local likelihood functions) are implicitly fused 
by the LC method in the process of calculating the JLF; they 
need not be explicitly transmitted between the sensors. There- 
fore, the LC method and our LC -based distributed (G)PFs are 
well suited for dynamic sensor networks. 

B. Simulation Setting 

In our simulations, the number of targets is P = 2 
unless stated otherwise. The system matrices G p and W p are 
identical for the two targets and given by lfT8l 



G p — 
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p = l,2. 



The variance of the driving noises u„ is given by a 2 = 
0.00035. Each of the two targets emits a sound of equal 
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(a) (b) (c) 



Fig. 1 . Example of a sensor network and communication topology, along with (a) a local likelihood function for one target, (b) a JLF for one target, and (c) 
a realization of the trajectories of two targets and the corresponding trajectories tracked by the LC-DPF. In (a), the square indicates the sensor for which the 
local likelihood is depicted. In (a) and (b), darker shading represents higher likelihood values and the cross indicates the position of the target. In (c), the 
stars indicate the start points of the target trajectories. 



amplitude A p = 10. The initial prior pdf /(x ) = 
A/"(Mo > Co) is different for the two targets, with p,^ = 
(36 36 -0.05 -0.05) T for target 1, p^ = (4 4 0.05 0.05) T 
for target 2, and C = diag{l, 1, 0.001, 0.001} for both 
targets. 

The network consists of K=25 acoustic amplitude sensors 
that are deployed on a jittered grid within a rectangular 
region of size 40 m x 40 m. Each sensor communicates with 
other sensors within a range of 18 m. The measurement noise 
variance is = 0.05 and the path loss exponent is k = 1. 

For LC, we approximate the measurement function 
h n ,k(x-n) in (l43l l by a polynomial (see d29l i) of degree R p = 2. 
This results in the following approximation of S n (z„,x n ) (cf. 
d32j): 

SVt( z m x n) 
25 4 

^E^Mtf'r (yi 1} r (4 2) r (^ r ■ 

fe=l r=0 

To obtain the approximation coefficients a n ,k,r needed for cal- 
culating the f5 n ,k,r{z n ,k) according to (FTJt and Pit , we use LS 
fitting as described in Section ITlI-CI The sums over all sensors 
in (l35l l are computed by average consensus algorithms using 
Metropolis weights J36). There are N c = ( 4 + 4 ) - 1 = 69 
consensus algorithms that are executed in parallel, each using 
1 = 8 iterations unless noted otherwise. The same remarks 
apply to the sums in fiTT i and (1421) . which are required by the 
R-LC-DGPF. The number of additional consensus algorithms 
employed by the R-LC-DGPF is N' c = 8 + 8 • 9/2 + 1 = 45. 

We compare the LC-DPF, LC-DGPF, R-LC-DGPF, CPF, 
and a centralized GPF (CGPF), which, similarly to the CPF, 
processes all sensor measurements at an FC. In addition, we 
consider the state-of-the-art consensus-based distributed PFs 
proposed (i) by Gu et al. in (211 (abbreviated GSHL-DPF), 
(ii) by Oreshkin and Coates in ED (OC-DPF), and (iii) by 
Farahmand et al. in |19i| (FRG-DPF). Unless stated otherwise, 
the number of particles at each sensor was J = 5000 for the 
LC-DPF, LC-DGPF, GSHL-DPF, and OC-DPF; J = 2000 for 
the FRG-DPF (this reduction is made possible by the adapted 



proposal distribution); and J' = 5000/25 = 200 for the R-LC- 
DGPF. The PF at the FC of the CPF and CGPF employed 
5000 particles. In the FRG-DPF [19), the rejection probability 
used for proposal adaptation was set to f3k = 0.02, and the 
oversampling factor was chosen as L = 10. 

As a performance measure, we use the n-dependent root- 
mean-square error of the targets' position estimate p n ,k, 
denoted RMSE„, which is computed as the square root of 
the average of \\p„\ — Pn > ' > || 2 over the two targets p = 1,2, 

all sensors k = 1.....25, and 5000 simulation runs. Here, 

(v) ~(v) 
pn denotes the position of target p and p n k denotes the 

corresponding estimate at sensor k. We also compute the 

average RMSE (ARMSE) as the square root of the average 

of RMSE^ over all 200 simulated time instants n. Finally, we 

assess the error variation across the sensors k by the standard 

deviation ctarmse of a fc-dependent error defined as the square 

root of the average of \\p„\ — pi^l) 2 ° ver the two targets 

p = 1,2, all 200 time instants n, and 5000 simulation runs. 

C. Simulation Results 

Fig- CD shows an example of a sensor network and commu- 
nication topology. For the case of a single target (P = 1), 
examples of the local likelihood function and of the JLF 
are visualized in Fig. [TJa) and (b), respectively. The local 
likelihood function is circularly symmetric because the mea- 
surement function /i n> fe(x n ) in (|43T > depends only on the 
distance between the target and the sensor. We can also see 
that the JLF is unimodal, which is an expected result since 
the JLF is the product of the local likelihood functions of 
all K = 25 sensors (see (Q}), all having circularly symmetric 
shapes as shown in Fig. [TJa) but different locations due to 
the different local measurements and the different distances 
between target and sensor (see (PPflO. Furthermore, we note that 
the nonlinearity of the local measurement functions /i„.fc(x n ) 
results in a non-Gaussian posterior (not shown in Fig. [TJ. 
For the case of two targets as described in Section IVIII-BI 
Fig. [Tic) shows a realization of the target trajectories and the 
corresponding tracked trajectories that were obtained at one 
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ARMSE [m] 


Track loss adjusted 
ARMSE [m] 


CT ARMSE [m] 


Track loss adjusted 
c ARMSE [m] 


Track loss 
percentage [%] 


Communication 
requirements 


LC-DPF 


0.6225 


0.5424 


0.0860 


0.0222 


0.95 


13800 


LC-DGPF 


0.6187 


0.5387 


0.0889 


0.0205 


0.7 


13800 


R-LC-DGPF 


0.5531 


0.5204 


0.0005 


0.0005 


0.46 


22800 


GSHL-DPF [2 11 


1.3022 


1.2841 


0.0032 


0.0032 


0.74 


8800 


OC-DPF [22.1 


0.9992 


0.8399 


0.0022 


0.0024 


1.1 


8800 


FRG-DPF 1 19| 


0.5553 


0.5335 








0.2 


400000 


CPF 


0.4975 


0.4975 









770 


CGPF 


0.5156 


0.5086 






0.18 


770 



TABLE I 

Estimation performance and communication requirements of the proposed consensus-based distributed PFs (LC-DPF, LC-DGPF, 

AND R-LC-DGPF), OF STATE-OF-THE-ART CONSENSUS-BASED DISTRIBUTED PFS (GSHL-DPF, OC-DPF, AND FRG-DPF), AND OF CENTRALIZED PFS 

(CPF AND CGPF). 



specific sensor by means of the LC-DPF. It can be seen that 
the target is tracked fairly well. Other sensors obtained similar 
results. 

Table U summarizes the estimation performance (ARMSE, 
track loss adjusted ARMSE, oarmse, track loss adjusted 
carmse, and track loss percentage) and the communication 
requirements of the proposed consensus-based distributed PFs 
(LC-DPF, LC-DGPF, and R-LC-DGPF), of the state-of-the-art 
consensus-based distributed PFs (GSHL-DPF, OC-DPF, and 
FRG-DPF), and of the centralized methods (CPF and CGPF). 
The "track loss percentage" is defined as the percentage of 
simulation runs during which the estimation error at time 
n = 200 exceeded 5m, which is half the average inter- 
sensor distance. Such simulation runs were excluded in the 
calculation of the "track loss adjusted" RMSE„, ARMSE, 
and cr armse- However, Table U presents also the ARMSE and 
carmse computed using all the simulation runs (including 
those with lost tracks). The "communication requirements" 
are defined as the total number of real numbers transmitted 
(over one hop between neighboring sensors) at one time instant 
within the entire network. For the centralized methods (CPF 
and CGPF), we used multi-hop routing of measurements and 
sensor locations from every sensor to the FC (located in one 
of the corners of the network). Furthermore, the estimates 
calculated at the FC are disseminated throughout the network, 
such that every sensor obtains the centralized estimate. 

It is seen from Table|I]that the track loss adjusted ARMSEs 
of the proposed distributed PFs are quite similar and that 
they are close to those of the centralized methods; they are 
slightly higher than that of FRG-DPF, slightly lower than that 
of OC-DPF, and about half that of GSHL-DPF. For FRG-DPF, 
carmse is zero, since max and min consensus algorithms are 
employed to ensure identical results at each sensor. Further- 
more, d armse is higher for LC-DPF and LC-DGPF than for 
R-LC-DGPF, GSHL-DPF, and OC-DPF. This is because R- 
LC-DGPF, GSHL-DPF, and OC-DPF employ a consensus step 
whereby Gaussian approximations of the partial/local posterior 
pdfs are combined to obtain a global posterior, thus achieving a 
tighter coupling between the sensors. By contrast, the local PFs 
of LC-DPF and LC-DGPF operate completely independently; 
only the JLF is computed in a distributed way using the 
LC scheme. Note, however, that the ARMSE and track loss 



adjusted ARMSE of LC-DPF and LC-DGPF are lower than for 
GSHL-DPF and OC-DPF. Finally, the track loss percentages of 
the proposed distributed PFs are below 1 % and similar to those 
of GSHL-DPF, OC-DPF, and FRG-DPF. As a consequence, 
the ARMSEs are generally very close to the track loss adjusted 
ARMSEs. 

The communication requirements of the distributed PFs are 
seen to be much higher than those of the centralized methods. 
This is due to our low-dimensional (scalar) measurements and 
the fact that each local likelihood function is parametrized only 
by the sensor location, i.e., three real numbers must be trans- 
mitted in one hop. For high-dimensional measurements and/or 
a different parametrization of the local likelihood functions, 
resulting in about 190 or more real numbers to be transmitted 
in one hop, the opposite will be true. Note that even when 
the consensus-based methods require more communications, 
they may be preferable over centralized methods because they 
are more robust (no possibility of FC failure), they require no 
routing protocols, and each sensor obtains an approximation 
of the global posterior (in the centralized schemes, each sensor 
obtains from the FC only the state estimate). It is furthermore 
seen that the communication requirements of the proposed 
distributed PFs are higher than those of GSHL-DPF and OC- 
DPF but much lower than those of FRG-DPF. Note, however, 
that the communication requirements of FRG-DPF depend on 
the number of particles and thus could be reduced by using 
fewer particles, whereas those of the other methods do not 
depend on the number of particles. (A setting with a lower 
number of particles will be considered later.) Finally, among 
the proposed distributed PFs, the communication requirements 
of R-LC-DGPF are higher by about 65% than those of LC- 
DPF and LC-DGPF. 

In Fig. |2 we compare the RMSE„ and track loss adjusted 
RMSE„ of the proposed LC-DGPF with that of CGPF and the 
state-of-the-art distributed PFs (GSHL-DPF, OC-DPF, FRG- 
DPF). In terms of track loss adjusted RMSE„ (Fig. Eft)), LC- 
DGPF outperforms GSHL-DPF and OC-DPF, and it performs 
almost as well as FRG-DPF and CGPF. The increase in 
RMSE„ over time in Fig. f2ja) is caused by the lost tracks. 

In Fig. [3] we compare the RMSE„ and track loss adjusted 
RMSE„ of LC-DPF (using eight consensus iterations) with 
that of CPF. As a performance benchmark, we also show 
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the results obtained by an impractical variant of LC-DPF in 
which the consensus algorithm is replaced by an exact, direct 
calculation of the sums in d35l l. The performance degradation 
of LC-DPF with exact sum calculation relative to CPF is 
caused by the LS approximation of the sensor measurement 
functions. The additional performance degradation of LC-DPF 
with eight consensus iterations relative to LC-DPF with exact 
sum calculation is due to the insufficiently converged consen- 
sus algorithms; it can be reduced by using more consensus 
iterations. In terms of the track loss adjusted RMSE„, both 
performance degradations are seen to be quite moderate. The 
track loss percentages were 0.95% for LC-DPF, 0.29% for 
LC-DPF with exact sum calculation, and 0% for CPF. 

Fig. [4] shows the track loss adjusted ARMSE of the 
proposed LC-DGPF and R-LC-DGPF versus the number / 
of consensus iterations. Here, R-LC-DGPF uses / consensus 
iterations in each one of its two consensus stages (i.e., / 
iterations to compute the sums in d35l l and / iterations each 
to compute the sums in (f4TT > and (l42l ). As a performance 
benchmark, the figure also shows the results for impractical 
variants of LC-DGPF and R-LC-DGPF using exact, direct 
calculation of the sums ( l35k (|4TT i. and ( l42l . It is seen that the 



performance of the impractical direct calculation is essentially 
achieved for / about 7 in the case of R-LC-DGPF and for 
/ about 10 in the case of LC-DGPF. Somewhat surprisingly, 
R-LC-DGPF outperforms LC-DGPF for up to 10 consensus 
iterations, i.e., the additional consensus algorithms used to cal- 
culate the sums in d4Tb and (l42l result in a better performance 
of R-LC-DGPF, in spite of the significantly reduced number 
of particles (200 instead of 5000). However, as the number 
of consensus iterations increases, both methods approach 
the performance of the respective "exact sum calculation" 
variant and LC-DGPF slightly outperforms R-LC-DGPF. This 
behavior can be explained as follows. The LC with a small 
number of consensus iterations is not completely converged, 
which means that the local information is not yet completely 
diffused throughout the network and the resulting approximate 
JLF does not yet contain the complete global information. 
The additional consensus stage of R-LC-DGPF then helps to 
further diffuse the local information. 

Finally, we consider a setting where each sensor in the 
distributed PF methods (LC-DPF, LC-DGPF, GSHL-DPF, OC- 
DPF, and FRG-DPF) as well as the FC in CPF and CGPF use 
only J = 400 particles, and consequently R-LC-DGPF uses 
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Fig. 4. Track loss adjusted ARMSE of the LC-DGPF and R-LC-DGPF versus 
the number / of consensus iterations, along with the track loss adjusted 
ARMSE of the impractical LC-DGPF and R-LC-DGPF variants with exact 
sum calculation. (R-LC-DGPF uses / consensus iterations for each sum 
calculation.) 



only J' = 400/25 = 16 particles per sensor. This reduction 
of the number of particles results in reduced communication 
requirements of FRG-DPF but not of the other methods as their 
communication requirements are independent of the number 
of particles. Table [II] summarizes the simulation results we 
obtained. A comparison with Table U shows that, as expected, 
the performance of all methods is degraded. Furthermore, the 
high ARMSE and track loss percentage values of LC-DPF, 
LC-DGPF, and OC-DPF can be viewed as signs of divergence. 
In the case of LC-DPF and LC-DGPF, high (Jarmse values 
indicate significant differences between the local particle rep- 
resentations of the global posterior; these differences reduce 
the effectiveness of the LS approximation in the LC scheme. 
In the case of OC-DPF, the divergence is due to the peaky 
functions (powers of local likelihoods functions) used in the 
weight update, which cause most of the particles to be located 
in regions of low likelihood. FRG-DPF performs well due to 
its use of adapted proposal distributions; its communication 
requirements are now closer to those of the other methods 
but still higher. R-LC-DGPF is seen to perform even slightly 
better with, at the same time, lower communication costs. As 
mentioned before, the additional consensus algorithms used 
by R-LC-DGPF lead to very similar particle representations 
of the local PFs across the network, with particles located in 
almost identical regions of the state space; this is evidenced by 
the low value of uarmse- Therefore, all sensors perform the LS 
approximation of their local likelihood functions in almost the 
same state space region, which moreover is the region where 
the particles of all sensors are located. Combining the local 
approximations using the LC scheme, we thus obtain a JLF 
approximation that is most accurate in that state space region. 
This explains the good tracking performance of R-LC-DGPF. 

IX. Conclusion 

For global estimation tasks in wireless sensor networks, 
the joint (all-sensors) likelihood function (JLF) plays a central 
role because it epitomizes the measurements of all sensors. We 



proposed a distributed, consensus-based method for computing 
the JLF. This "likelihood consensus" method uses iterative 
consensus algorithms to compute, at each sensor, an approxi- 
mation of the JLF as a function of the state to be estimated. 
Our method is applicable if the local likelihood functions of 
the various sensors (viewed as conditional probability density 
functions of the local measurements) belong to the exponential 
family of distributions. This includes the case of additive Gaus- 
sian measurement noises. The employed consensus algorithms 
require only local communications between neighboring sen- 
sors and operate without complex routing protocols. 

We demonstrated the use of the likelihood consensus 
method for distributed particle filtering and distributed Gaus- 
sian particle filtering. At each sensor, a local particle filter 
computes a global state estimate that reflects the measure- 
ments of all sensors. The approximate JLF provided by the 
likelihood consensus method is used for updating the parti- 
cle weights of each local particle filter. A second stage of 
consensus algorithms can be employed to significantly reduce 
the complexity of the distributed Gaussian particle filter. We 
applied the proposed distributed particle filters to a multiple 
target tracking problem and demonstrated experimentally that 
their performance is close to that of the centralized particle 
filters. Compared to three state-of-the-art distributed particle 
filtering schemes, our methods typically achieve a comparable 
or better estimation performance, while the communication 
requirements are somewhat higher in two cases and much 
lower in one case. 

We finally note that the proposed distributed Gaussian 
particle filter can be extended to a consensus-based, distributed 
implementation of the Gaussian sum particle filter proposed 
in ll46l . Furthermore, an extension of the likelihood consensus 
method to general local likelihood functions (i.e., not neces- 
sarily belonging to the exponential family) has been presented 

in ea. 
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